For  conversion  of  SI  metric  units  to  U.S./British 
customary  units  of  measurement  consult  ASTM 
Standard  1380,  Metric  Practice  Guide,  published 
by  the  American  Society  for  Testing  and 
Materials,  1916  Race  St.,  Philadelphia,  Pa.  19103. 


Cover  Accretion  of  rime  ice,  formed  from  sub 
freezing  cloud  droplets,  on  a suspender 
cable  on  the  summit  of  Ml  Washington. 
New  Hampshire,  in  April  1978  The  cable 
was  instrumented  to  measure  the  addi 
tional  loading  placed  on  it  by  wind  and 
gravity  forces  acting  on  the  attached  ice 
(Photograph  by  Ion  L ingel,  Ml  Wash- 
ington Observatory.) 


I 


CRREL  Report  79-4 


Computer  modeling  of  atmospheric 
ice  accretion 


S.F.  Ackley  and  M.K.  Templeton 


Prepared  for 

DIRECTORATE  OF  MILITARY  PROGRAMS 
OFFICE,  CHIEF  OF  ENGINEERS 
Bv 

UNITED  STATES  ARMY 
CORPS  OF  ENGINEERS 

COLD  REGIONS  RESEARCH  AND  ENCINEERINC  LABORATORY 
HANOVER,  NEW  HAMPSHIRE,  U S A. 


Approved  for  pubJu  release,  distribution  unlimited 


SECURITY  CLASSIFICATION  OF  THU  PACE  (When  Dele  entered) 

J REPORT  DOCUMENTATION  PAGE 


12.  GOVT  ACCESSION  NO. 


CRREL 


l».  TITLE  (and  Subtitle) 


I j COMPUTER  MODELING  OF  ATMOSPHERIC  JCE  ACCRETION, 


READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

3.  RECIPIENT'S  CATALOG  NUMBER 


S.  TYPE  OF  REPORT  A PERIOD  COVERED 


S.  PERFORMING  ORG.  REPORT  NUMBER 


17.  AUTHOR!’.; 


8.  CONTRACT  OR  GRANT  NUMBERfa.) 


mt  F ^Ackley  aMtM.K^Yempleton 


».  PERFORMING  ORGANIZATION  NAME  AND  ADORESS 

U.S.  Army  Cold  Regions  Research  and  Engineering  Laborato 
Hanover,  New  Hampshire  03755  S' 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Directorate  of  Military  Programs 
Office,  Chief  of  Engineers 
Washington,  D.C.  20314 

14.  MONITORING  AGENCY  NAME  a ADORESSffr  dlttarant  I 


i Controlling  Olllea) 


10.  PROCmU-BCEMENT.  PROJECT,  TASK 
AREA  nMiV 

DA  Proiect|4A161 1J32AT24  j 
DA  Project  4A7670AT42  ~ 

report  date 

f V — 7 

- / Mare>M>79  / 

NUMBER  OF  PAGES  i 

39 r. 

IS.  SECURITY  CLASS.  ( ol  thla  report) 

Unclassified 

ISa.  DECLASSIFICATION/DOWNGRADING 
SCHEDULE 


I IS.  DISTRIBUTION  STATEMENT  (ol  thla  Raport) 


Approved  for  public  release;  distribution  unlimited. 


I 17.  DISTRIBUTION  STATEMENT  (ot  the  abatract  antarad  In  Block  30,  II  dlttarant  from  Raport ) 


IB.  SUPPLEMENTARY  NOTES 


I IS.  KEY  WORDS  (Continue  on  rararaa  at  da  It  nacaaaary  mtd  Identity  by  block  number ) 


Cold  regions  Ice 

Computerized  simulation  |ce  formation 
C°ntr°l  Mathematical  models 

Deicing  systems  Test  methods 

Drops Clouds 

ABSTRACT  (CeeUmae  amteeerwa  atda  m meoaoeay  mag  Identity  by  block  number} 

A computer  model  is  described  to  compute  the  amount  of  ice  accretion  on  an  object  under  a variety  of  initial  condi- 
tions. Numerical  techniques  are  best  applied  to  these  problems  because  of  time  dependent  effects  governing  the 
amount  of  ice  collected  and  also  the  variety  of  initial  conditions  that  can  lead  to  ice  accumulation.  The  helicopter 
rotor  icing  problem  adds  an  additional  complexity  since  the  velocity  along  the  rotor  blade  varies  over  a wide  range 
strongly  affecting  the  amounts  of  ice  collected  at  different  blade  positions.  The  physics  of  ice  accretion  is  reviewed 
and  the  accounting  for  time-dependence  in  the  computer  model  is  described.  Some  model  results  are  presented  and 
indicate  the  dependence  of  ice  accretion  on  velocity,  droplet  sizes,  cloud  liquid  water  content,  and  temperature  for  a 
cylindrical  object  of  constant  size.  


DD/jSTyyUn 


KDfTION  or  f NOVI 


037  3-00 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  ( 


PREFACE 

This  report  was  prepared  by  S.F.  Ackley,  Research  Physicist,  Snow  and  Ice  Branch, 
Research  Division  and  M.K.  Templeton,  Math  Aide,  Snow  and  Ice  Branch  (now  at 
California  Institute  of  Technology). 

The  study  was  conduct'  j under  DA  Project  4A1 61 1 02AT24,  Adhesion  and  Physics 
of  Ice  and  DA  Project  4 A / ''7  3AT42,  Prevention  and  Control  of  Ice  Adhesion. 

Mr.  Walter  B.  Tucker  III  a id  Mr.  L.  David  Minsk  reviewed  the  technical  content  of 
this  report. 

The  authors  wish  to  thank  Dr.  G.  Ashton,  Chief,  Snow  and  Ice  Branch  and  Program 
Manager  for  CRREL’s  Ice  Adhesion  Program  for  his  continued  support  and  encourage- 
ment of  this  work. 


CONTENTS 


Page 


Abstract j 

Preface ii 

Introduction 1 

Ice  accretion  parameters 2 

Physics  of  ice  accretion 2 

1 . Interaction  between  water  droplets  and  flow  field 2 

2.  Time  dependence  in  droplet  trajectories 4 

3.  Thermodynamic  processes  at  the  surface  of  the  object 4 

4.  Time  dependence  in  the  thermodynamics 6 

5.  Time  dependence  of  lesser  order 6 

Numerical  ice  accretion  model 7 

1.  Major  subroutines 7 

2.  Options  for  droplet  size  variations 1 1 

3.  Option  for  the  helicopter  rotor  case 1 5 

Results 16 

Conclusions  and  future  studies 20 

Literature  cited 21 

Appendix  A:  Computer  program 23 


ILLUSTRATIONS 

Figure 

1 . Velocity  potential  </>  and  stream  function  ^ around  a circular  cylinder  of  radius  a em- 

bedded in  a fluid  with  free  stream  velocity  u 3 

2.  Air  streamlines  and  droplet  trajectories  with  respect  to  a right  circular  cylinder 4 

3.  Schematic  indicating  the  five  terms  used  in  the  heat  balance  and  whether  they  are 

carrying  heat  toward  or  away  from  the  freezing  surface 6 

4.  The  fraction  of  accreted  water  lost  for  various  temperatures  and  average  radii  of  accret- 

ing objects 7 

5.  Block  diagram  of  the  numerical  model  for  calculating  ice  accretion  8 

6.  Initial  variables  input  in  the  run  control  sequence  of  the  numerical  model 8 

7.  Subroutines  Profile  and  Ice  Oen  used  to  calculate  the  changes  in  object  profile  and  ice 

density  from  the  output  of  the  icing  thermodynamics 9 

8.  Subroutines  controlling  computation  of  the  droplet  trajectory,  collection  efficiency  and 

trajectory  control,  and  free  stream  velocity 1 0 

9.  Based  on  the  ambient  conditions,  the  calculation  of  mass  accreted  is  completed  by  in- 

corporating the  thermodynamics  at  the  surface 1 3 

1 0.  Forms  of  the  Erlang  distribution  for  various  values  of  the  parameter  k 14 

1 1 . (T op)  The  change  in  profile  dimension  at  50-s  intervals  by  ice  accretion 1 5 

(Bottom)  The  surface  temperature  of  the  front  half  cylinder  as  a function  of  time 1 5 

1 2.  Numerical  simulation  of  helicopter  rotor  blade  icing  showing  the  surface  temperature, 

leading  edge  ice  thickness,  collection  efficiency  and  fraction  accreted  as  a function  of 
velocity  for  the  initial  conditions  shown. 16 

1 3.  Droplet  trajectories  and  ice  profile  changes  for  a Gaussian  droplet  size  distribution 1 7 

14.  Accreted  ice  thickness  and  temperature  for  the  ambient  conditions  indicated 18 

15.  Comparison  between  experimental  data  and  model  simulation  of  the  same  conditions 19 

16.  Magnitude  of  the  individual  thermodynamic  terms  as  a function  of  velocity  for  the  simu- 

lation shown  in  Figure  1 2 20 


iii 


I 


i 


I 


COMPUTER  MODELING  OF 
ATMOSPHERIC  ICE  ACCRETION 


S.F.  Ackley  and  M.K.  Templeton 

INTRODUCTION 

Work  on  ice  accretion  from  the  atmosphere  has 
followed  two  distinct  lines,  hailstone  growth  and  air- 
craft icing.  Most  recently,  investigations  of  hailstone 
growth  have  provided  the  justification  for  a general 
understanding  of  ice  accretion  processes.  The  concept 
is  that  damaging  hail  grows  in  several  regions  within 
convective  cumulus  clouds.  Through  examination  of 
the  structure  of  hailstones,  information  on  the  ambient 
conditions  within  the  clouds  can  be  obtained,  specif- 
ically liquid  water  content  (Iwc),  droplet  radius  (f?drop) 
and  air  temperature  (T0).  The  hailstones  can  then  be 
used  to  interpret  and  model  cloud  structures,  and 
possibly  to  provide  information  on  structure  modifica- 
tion by  cloud-seeding  methods  for  preventing  and 
slowing  down  hail  growth.  Macklin  et  al.  (1976,  1977), 
List  et  al.  (1971),  and  Knight  and  Knight  (1968)  pro- 
vided the  conceptual  background  for  this  application 
and  described  several  recent  efforts  to  obtain  this  in- 
formation. 

Early  ice  accretion  work  (e.g.  Langmuir  and  Blodgett 
1946)  emphasized  its  application  to  airplane  icing.  Icing 
loads  affected  aircraft  performance,  making  it  necessary 
to  provide  ice  detection  and  deicing  systems  (Tribus 
1943,  Messinger  1953).  As  aircraft  evolved  to  types 
that  flew  higher  and  faster,  however,  research  in  this 
area  declined,  since  major  icing  problems  occur  at 
altitudes  below  20,000  ft  (6100  m)  and  at  air  speeds 
below  300  mph  (134  m/s),  primarily  in  clouds  ranging 
from  about  0 to  -20°C  in  temperature.  The  faster  and 
higher  flying  commercial  and  military  planes  of  today 
pass  quickly  through  these  conditions.  Jet  aircraft  also 
have  power  available  to  overcome  icing  on  engine  in- 
takes and  other  critical  areas  by  heating,  so  recent 
efforts  have  been  directed  toward  finding  engineering 
solutions  to  typical  or  worst  cases  rather  than  toward 
continued  research  into  the  physics  of  the  ice  accretion 
process. 

The  fundamental  physics  indicated  that  particular 
cases  could  be  formulated,  but  that  the  computations 
necessary  to  give  a complete  analysis  were  complex. 

For  example,  cumbersome  analog  devices  requiring 
several  persons  working  simultaneously  were  necessary 


to  obtain  the  initial  water  droplet  trajectories  of 
impaction.  These  initial  trajectories  were,  of  course, 
only  one  aspect  of  the  total  problem  (Brun  1957, 
Langmuir  and  Blodgett  1946). 

Research  on  helicopters  and  the  possibility  of  making 
them  capable  of  operation  in  any  weather,  including 
icing  conditions,  has  proceeded  sporadically  but  with 
continued  interest  since  helicopter  flight  is  limited 
mainly  to  the  lower  altitudes  and  air  speeds  where  icing 
is  more  frequently  encountered.  Stallabrass  (1957, 

1958)  conducted  experiments  on  helicopters  hovering 
under  various  icing  conditions,  using  a water  spray 
stand  to  simulate  the  icing  cloud.  Using  the  formulation 
developed  from  earlier  aircraft  work,  he  found  that  the 
icing  accretion  rate  on  the  main  rotor  blades  could  be 
simulated  by  applying  the  fixed  velocity  calculation 
developed  by  Messinger  (1953)  for  fixed-wing  aircraft 
to  each  segment  along  the  rotor  blade  characterized  by 
its  velocity: 

vi  = (1) 

where  w is  the  angular  velocity  of  the  rotor  blade  and 
Zj  is  the  distance  from  the  center  of  the  hub  (axis  of 
rotation)  to  the  segment  of  interest.  Good  agreement 
was  found  between  the  experimentally  obtained  ice 
shapes  and  those  predicted  using  the  calculations  from 
the  earlier  work  on  aircraft  at  a single  velocity  (Dickey 
1952). 

The  aim  of  this  report  is  to  demonstrate  a method 
for  calculating  the  quantity  of  ice  that  will  accrete  on 
fixed-wing  aircraft,  helicopters,  ships,  power  lines, 
radomes,  towers  and  other  structures  under  a given  set 
of  conditions. 

Although  much  of  the  previous  work  has  been 
fundamentally  correct  in  formulating  the  time-dependent 
process  of  ice  accretion,  one  of  the  problems  has  been 
the  exclusion  of  time  dependence  in  the  application  of 
the  formalism.  In  order  to  estimate  ice  accretion  mass 
confidently  over  a wide  range  of  conditions,  the  in- 
fluence of  this  time  dependence  must  be  accounted 
for. 

In  this  paper,  we  reexamine  the  fundamentals  of  the 
ice  accretion  process,  and  point  out  those  features  in 


the  governing  relationships  for  ice  accretion  that  require 
special  treatment  because  they  implicitly  contain  time 
dependence  leading  to  feedback  effects.  We  then  show 
how  the  time-dependent  formalism  is  solved  using 
numerical  methods  to  more  correctly  account  for  these 
feedback  effects,  and  include  them  in  calculations  of 
ice  accretion  on  objects,  given  the  atmospheric  con- 
ditions and  object  parameters.  The  computer  model 
is  useful  as  a simulation  tool  to  provide  input  for  de- 
sign purposes,  such  as  aircraft  and  helicopter  deicing 
design,  or  ground-based  structural  applications,  such 
as  power  line  or  tower  icing  prediction,  or  as  a research 
tool  to  examine  the  sensitivity  of  the  icing  process  to 
variations  in  the  input  variables,  either  singly  or  in 
concert.  These  simulations  may  assist  in  the  under- 
standing of  limited  laboratory  and  field  testing,  lend- 
ing some  generality  to  the  ice  accretion  processes 
observed  there,  and  extending  this  information  to 
applications  where  data  are  limited,  for  example, 
power  line  icing  in  remote  locations,  and  helicopter 
icing  in  the  extreme  cases  that  are  currently  untested. 


ICE  ACCRETION  PARAMETERS 

Atmospheric  ice  accretion  is  caused  principally  by 
the  freezing  of  supercooled  water  droplets.  The 
droplets  may  freeze  after  they  are  carried  to  a collision 
with  an  object  by  airflow  (ground-based  structure),  or 
as  an  object,  such  as  an  airplane  or  helicopter,  passes 
through  them.  Supercooled  water  is  necessary  because 
ice  particles  alone  do  not  stick  on  impact  and,  there- 
fore, do  not  cause  an  ice  accretion  problem.  Recently 
there  have  been  some  indications  that  mixed  clouds, 
containing  both  ice  particles  and  supercooled  water 
droplets,  can  cause  a serious  ice  accretion  problem, 
especially  for  helicopters  (Stallabrass,  pers.  comm.). 

This  problem  ostensibly  occurs  because  the  ice  particles, 
in  the  presence  of  supercooled  water,  adhere  to  the 
surface  of  the  object,  and,  because  they  are  already 
frozen,  do  not  require  the  heat  transfer  that  an  un- 
frozen droplet  of  similar  size  requires  to  compensate 
for  the  large  latent  heat  of  freezing.  However,  icing 
tunnel  tests  have  indicated  that  the  amount  of  stick- 
ing by  ice  crystals  is  lower  than  expected,  indicating 
that  mixed  conditions  may  not  be  as  critical  as  was 
first  assumed. 

In  the  first  case,  we  will  consider  only  water  droplets, 
since  the  effects  of  mixed  conditions  have  not  been 
studied  sufficiently  to  allow  a quantitative  formulation. 
The  primary  sources  of  icing  conditions  are  clouds  and 
fogs  since  the  smaller  droplets  present  in  them  can 
sustain  substantial  supercooling  over  a wide  temperature 


range.  Raindrops  are  larger  and  therefore  cannot 
undergo  the  same  degree  of  supercooling  (Mason  1971 ). 
However,  the  formulation  can  be  applied  to  clouds,  fog, 
rain,  sea  spray,  and  spray  of  man-made  origin,  e.g.  power 
plant  cooling  towjr  fog,  as  long  as  the  correct  parameters 
can  be  specified.  The  six  parameters  necessary  for 
quantifying  the  problem  can  be  conveniently  divided 
into  three  properties  of  the  atmosphere-the  ambient 
temperature  T0,  the  liquid  water  content  Iwc,  and  the 
droplet  radius  distribution  /?drop  -and  three  of  the 
accreting  object-its  cross-sectional  diameter  (2a), 
velocity  u and  shape,  e.g.  cylinder,  airfoil,  etc.  The 
interaction  of  these  parameters  will  be  identified  in  the 
next  section. 


PHYSICS  OF  ICE  ACCRETION 

If  we  move  an  object  through  a supercooled  cloud, 
the  rate  of  ice  accretion  depends  on  two  factors.  First, 
it  depends  on  the  kinematic  interaction  between  the 
object  and  the  droplets  in  the  cloud.  This  relation 
determines  whether  particular  droplet  sizes  are  captured. 
Second,  it  depends  on  the  thermodynamic  processes 
at  the  surface  of  the  object.  Here,  a balance  is  obtained 
between  the  rate  of  heat  release  necessary  to  freeze  all 
or  part  of  the  impinging  water  and  the  rate  at  which 
heat  can  be  carried  away  from  the  surface  into  the  flow 
field  by  the  processes  of  convection,  evaporation  and 
radiation.  We  first  consider  the  kinematic  interaction. 

1.  Interaction  between  water  droplets 
and  flow  field 

Normally,  the  speeds  at  which  icing  occurs  are  such 
that  the  air  around  the  object  can  be  treated  as  an  ideal, 
incompressible  fluid.  The  velocity  field  exhibits  con- 
tinuity, and  mathematically  (Sokolnikoff  and  Redheffer 
1966)  satisfies  the  condition  for  the  application  of 
potential  theory  (i.e.  Vxu  = 0,  V • u = 0).  The 
velocity  can  therefore  be  solved  for  at  any  point  around 
the  object  by  solution  of  a complex  potential  of  the 
form: 

F = <t>  + i*  (2) 

where  0 is  the  velocity  potential  and  0 the  stream 
function,  and  the  curves  of  0,  0 are  orthogonal  func- 
tions defining  the  complex  flow  field.  The  velocity  of 
the  fluid  at  any  point  is  given  by 

u0  = V0.  (3) 
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Because  of  the  orthogonality  between  <p  and  i//,  the 
velocity  is  tangent  to  the  curves  ^ = constant.  These 
are  the  paths  or  streamlines  of  the  fluid  elements.  We 
consider  a circular  cylinder  as  an  example.  Figure  1 
shows  the  streamlines  around  such  a circular  cylinder, 
with  flow  perpendicular  to  the  long  axis.  Two  con- 
sequences of  complex  potential  flow  theory  are  that 
the  boundary  of  the  object  coincides  with  one  of  the 
streamlines  (\l/  = const,  at  the  boundary),  and  the 
flow  satisfies  Laplace’s  equation,  V2/7  = 0,  in  the 
entire  plane.  For  objects  of  simple  cross-sectional 
geometry,  such  as  circles  or  ellipses,  the  dimensionless 
velocity  at  any  point  can  be  determined  straightfor- 
wardly if  we  know  the  free  stream  velocity  u (i.e.  at 
great  distances  from  the  object)  and  the  characteristic 
dimensions  of  the  object— for  the  cylinder,  its  radius  a 
(Milne-Thomson  1 960). 

To  determine  the  interaction  between  the  droplets 
and  the  flow  field,  we  integrate  the  equation  of  motion 
for  the  droplet  in  the  flow  field.  Following  Brun  (1957), 
the  assumptions  for  this  formulation  are:  a)  the  stream- 
lines determined  for  clear  air  are  valid,  i.e.  there  are  not 
enough  droplets  to  perturb  the  flow,  b)  the  gravity 
forces  are  much  less  than  the  inertial  forces  and  may 
be  neglected,  c)  the  pressure  forces  on  a droplet  are 
equivalent  to  those  on  an  equal  volume  of  air  at  the 
same  location,  and  consequently  may  also  be  neglected 
because  water  has  a much  greater  density  than  air. 

The  motion  of  the  droplet,  therefore,  results  pri- 
marily from  the  inertial  and  viscous  drag  forces  as  it 
drifts  with  the  fluid  streamline.  The  drag  force  is 
proportional  to  the  velocity  and  the  drag  coefficient 
cd  that  the  droplet  sees  at  any  given  time.  Newton’s 
second  law  in  dimensionless  form  for  the  droplet  (after 
Langmuir  and  Blodgett  1 946)  therefore  becomes,  for 
the  circular  cylinder  case: 


Figure  1.  Velocity  potential  <p  and  stream  function  <1/ 
around  a circular  cylinder  of  radius  a embedded  in  a 
fluid  with  free  stream  velocity  u. 


Fiere  ud  is  the  dimensionless  velocity  of  the  droplet  at 
any  point,  u0  the  dimensionless  streamline  velocity  at 
the  same  point,  Re  the  Reynolds  number  = 2/?drop 
Pair  l (^d  “ w0)o]  Ip,  cd  the  drag  coefficient,  and  k the 
inertial  parameter  (analogous  to  mass  in  dimensionless 
form)  = 2pw/?dr0(,  (<y/9po). 

The  drag  coefficient  cd  of  a spherical  drop  is  also  a 
function  of  Reynolds  number,  and  its  dependence  is 
given  by  Beard  and  Pruppacher  (1969)  as: 

cd  = 1 +0.102  Re-955  0.2<Re<2  ] 


cd  = 1+0.115  Re-802 


2<Re<21 


- (5) 


cd  = 1 +0.189  Re-632  21  < Re< 200 


The  calculation  of  droplet  trajectories  involves  con- 
siderable computation  time,  since  the  equation  of 
motion  has  to  be  integrated  numerically  with  respect 
to  time  to  fully  assess  the  degree  to  which  the  droplet 
deviates  from  the  fluid  streamlines.  Langmuir  and 
Blodgett  (1946)  compiled  these  computations  into 
curves  of  the  collection  efficiency  E as  a function  of 
two  dimensionless  parameters,  each  depending  on  the 
free  stream  velocity,  droplet  radius  and  cylinder  radius. 
Collection  efficiency  is  defined  as  the  distance  from 
the  centerline  of  the  cylinder  axis  divided  by  the 
cylinder  radius  that  a droplet  can  be  and  still  be 
captured  by  the  cylinder,  i.e.  the  far  field  position  of 
the  droplet  trajectory  that  is  tangent  to  the  cylinder 
(Fig.  2).  It  also  refers  to  that  fraction  of  the  total 
possible  droplets  in  the  path  of  the  object  that  collected. 
For  example,  a collection  efficiency  of  0.5  indicates  that 
half  the  droplets  of  a given  size  would  be  collected  in 
the  cross-sectional  area  swept  out  by  the  cylinder  in 
the  flow.  Linder  natural  conditions  there  is  usually  a 
variety  of  droplet  sizes,  so  the  computations  become 
lengthy  since  an  equation  of  motion  has  to  be  integrated 
separately  for  each  size  category.  The  total  collection 
efficiency  is  then  the  sum  of  the  collection  efficiencies 
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Figure  2.  Air  streamlines  and  droplet  trajectories  with  respect  to  a right  circular  cylinder.  Collection  ef- 
ficiency E is  the  ratio  of  the  initial  far-field  y location  y0  of  the  droplet  trajectory  that  is  just  tangent  to  the 
radius  of  the  cylinder  a. 


for  each  size  times  the  fraction  of  the  total  Iwc  repre- 
sented by  that  size. 

2.  Time  dependence  in  droplet  trajectories 

The  first  of  the  time-dependent  effects  that  have 
not  been  treated  in  detail  in  earlier  work  is  that  of 
droplet  trajectory.  As  ice  accretes  on  the  front  surface 
of  the  cylinder  it  is  apparent  that  the  cylinder  cannot 
maintain  its  initial  circular  shape.  Therefore,  the  stream- 
line characterized  by  the  boundary  of  the  object  is  no 
longer  the  same  as  it  was  initially,  and  the  velocity 
field  in  other  regions  is  proportionally  distorted  to 
maintain  continuity.  From  eq  5,  the  change  in  velocity 
of  the  droplet  depends  on  the  streamline  velocity  in 
that  location.  The  trajectory  and  velocity  of  a droplet 
of  a given  size  will  therefore  change  with  time  from 
their  values  prior  to  ice  accretion,  thereby  changing 
the  collection  efficiency  E as  ice  accretes.  This  effect 
has  not  been  included  systematically  before  because 
of  the  length  of  the  computations  required  to  get  the 
initial  pre-icing  values  of  the  collection  efficiency. 

With  a digital  computer  of  sufficient  power,  however, 
this  effect  can  be  taken  into  consideration  as  we  show 
in  a later  section. 

3.  Thermodynamic  processes  at  the  surface 
of  the  object 

In  general,  the  mass  rate  of  water  arriving  at  the 
surface  of  the  object  (per  unit  length)  can  be  written 
as: 


where  E is  the  collection  efficiency,  dependent  on  the 
velocity,  the  droplet  size,  and  the  object’s  cross- 
sectional  diameter  and  characteristic  dimensions;  iwc 
is  the  volumetric  liquid  water  content  (kg/m3);w  is 
the  free  stream  velocity  (m/s);  and  2 a is  the  object’s 
cross-sectional  diameter  normal  to  the  flow  (m).  How- 
ever, the  mass  of  ice  accreted  will  not  in  general  be 
equal  to  this  quantity  if  heat  transfer  processes  cannot 
adequately  freeze  all  of  the  accreted  water. 

In  this  formulation  it  is  assumed  that  1)  radiation 
terms  are  excluded  since  they  are  small  relative  to 
convective  and  evaporative  heat  flux,  2)  near  the  ob- 
ject, the  droplets  have  the  same  velocity  as  air  at  that 
location  (»«0),  3)  the  droplets  are  at  thermodynamic 
equilibrium  at  temperature  T0,  the  ambient  air 
temperature,  and  4)  the  boundary  layer  thickness  is 
small  so  the  droplet  passes  through  the  boundary  layer 
and  strikes  the  surface  with  the  temperature  T0  and 
velocity  u that  it  had  in  the  flow.  Recent  work  (List 
1 977,  List  et  al . 1 976, ) oe  et  al.  1 976)  has  indicated 
some  problems  with  these  assumptions,  especially  with 
the  boundary  layer  interaction  and  the  droplet’s  thermal 
state  at  high  liquid  water  contents  (>0.002  kg/m"*). 

But,  as  we  indicate  later,  inclusion  of  different  thermo- 
dynamics and  assumptions  can  be  accommodated  in 
numerical  solutions.  We  employ  those  just  described 
because  of  their  long  usage  and  their  verification,  in 
general  terms,  by  experimental  observations. 

The  heat  balance  of  the  surface,  controlled  by  the 
properties  of  the  flow,  is  taken  as  the  sum  of  the 
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convective  and  evaporative  fluxes,  and  the  aerodynamic 
heating  of  the  incoming  air: 

<?surf  = -4(r-r0)-^-+^(f'-f0)|-  P) 
^•pair  ^Lpair 


New  terms  are  H,  the  convective  heat  transfer  coef- 
ficient (defined  below),  R , the  surface  recovery  factor, 
cp  ajr,  the  heat  capacity  of  air,  L0,  latent  heat  of  water 
vaporization,  E' , vapor  pressure  of  water  at  surface 
temperature  T,  E0,  vapor  pressure  of  water  at  the 
ambient  temperature  T0,  an  dp,  atmospheric  pressure. 
The  factor  0.6/p  converts  the  mass  flux  from  weights 
of  water  vapor  to  the  more  convenient  vapor  pressure 
form,  and  uses  the  relationship  of  latent  and  sensible 
heat  transfer  coefficients  of  Hm/H  s [F(Le)/cp  air] 
where  F(Le)  (a  function  of  the  Lewis  number  Le  (ratio 
of  the  diffusion  coefficient  to  thermal  diffusivity ) ] is 
approximately  1 . Therefore,  Hm  a (H/cp  ajr),  as  indi- 
cated in  eq  7. 

The  heat  transfer  coefficient  per  unit  area,  H0,  is 
taken  from  Bosch’s  formulation  (Boelter  et  al.  1965): 


H0  = 31.01 


tv-J0 

(2j)-44 


(P:  2?3J6\-S6 
\ 1014 F*  ) 


Here  F*  is  the  film  temperature  = (T0  + 7")/2.  Initially 
the  surface  temperature  T is  taken  as  the  recovery  tem- 
perature T = T0  + (Rul/c p air)  where  Ru$/cp  air  is 
the  aerodynamic  heating  effect  (in  °C)  from  com- 
pression of  the  non-icing  flow  at  the  stagnation  point. 
The  heat  transfer  coefficient  H is 


H = H0A 


where  A is  the  surface  area  over  which  the  identified 
processes  are  taking  place,  taken  as  the  front  half 
cylinder  for  our  example. 

The  heat  given  to  the  blade  by  the  incoming  drop- 
lets (mass/n  per  unit  time)  is  dependent  on  whether 
none,  some,  or  all  of  the  water  is  accreted  as  ice,  that 
is,  the  three  cases  of  surface  temperature:  1)  T> 0°C, 

2)  T = 0°C,  3)  T< 0°C. 

Case  1.  r>0°C 

In  this  case,  the  heat  flux  is  due  to  the  temperature 
difference  of  the  collected  water  plus  the  conversion 
of  the  droplet’s  kinetic  energy  into  heat: 

<H" 


Case  2.  T=0°C 

In  this  case,  a fraction  Fof  the  accreted  water  is 
converted  to  ice,  adding  the  latent  heat  of  fusion  L to 
the  heat  flux: 

<?2  = ™ph2o  • (11) 

2 cp  h2o  2cp  h2o 


Case  3.  7-  0°  C 

Here  al1  the  accreted  water  is  frozen  as  ice,  i.e.  F = 

1 and  the  first  term  in  the  following  equation  accounts 
for  the  difference  in  specific  heat  between  the  water 
and  ice  at  temperature  T: 


<?3  = me. 


h2o  (T^-To-rL--2T—\{u) 

2 \ CPh2°  cP  H20  2cpH20 J 


The  thermodynamic  steady-state  balance  requires 
that  the  heat  taken  by  the  flow,  C?surf,  is  balanced  by 
the  appropriate  value  (?,,  Q2  or  Q3.  In  order  to  choose 
between  these  cases,  Brun  (1957)  and  Messinger  (1953) 
have  taken  the  heat  balance  for  Case  2 and  solved  for 
the  frozen  fraction  F.  If  F is  between  the  limits  0 and 
1 then  Case  2 is  the  appropriate  one.  If  F>  1 then  Case 
3 is  applicable,  and  if  F< 0 then  Case  1 applies. 

In  a similar  manner,  the  heat  balance  as  shown 
schematically  in  Figure  3 (after  Messinger  1953)  can 
be  taken  for  any  case,  and  the  residual  heat  in  a short 
time  interval  At 

Q = «?surf-<?i)Af  '=  I-2,3  (1 

can  then  be  used  to  modify  the  surface  temperature  by 

■r  _ 7cp  m -blade  + Q fi 

.iaw  “ , • ' 

cp  m-blade 


■'p  m-blade  cp  m-blade 


+ m-F-c„ 


i.e.  the  total  heat  capacity  of  the  object  (blade)  plus 
the  heat  capacity  of  any  accreted  ice.  The  new  tem- 
perature 7"new  can  then  be  checked  to  see  if  it  is 
above  or  below  0°C  and  the  appropriate  thermodynamic 
branch  ((?,,  <?2,  (?3)  selected.  In  going  from  above  or 
below  0°C,  the  option  T = 0°C  is  selected  and  the 
appropriate  balance  chosen  based  on  the  value  of  the 
accreted  fraction  F.  F must  lie  within  the  range  0 <F<  1 
in  order  for  the  T - 0 thermodynamics  to  be  chosen. 
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Figure  3.  Schematic  indicating  the  five  terms 
used  in  the  heat  balance  and  whether  they  are 
carrying  heat  toward  ( +)  or  away  (-)  from  the 
freezing  surface.  Convective  flux  usually  takes 
heat  away  from  the  surface,  but  at  high  veloc- 
ities adiabatic  compression  of  the  flow  domi- 
nates and  the  convective  heat  flux  is  then 
positive  toward  the  surface. 


4.  Time  dependence  in  the  thermodynamics 

As  with  the  droplet  trajectories,  the  time  dependency 
of  the  thermodynamics  associated  with  ice  accretion 
basically  results  from  the  distortion  of  the  object’s  size 
and  shape  as  the  ice  accretes. 

Equation  6 for  the  mass  rate  of  water  collected  con- 
tains this  dependence  in  the  collection  efficiency  term, 
which  changes  as  the  flow  field  responds  to  the  in- 
creasing accretion.  The  heat  flux  from  the  surface  into 
the  flow  (eq  7)  varies  with  the  change  in  heat  transfer 
coefficient  (eq  9),  which  in  turn  varies  through  the 
changes  in  surface  area  and  shape  with  time.  Since  the 
heat  fluxes  Q,  additionally  depend  on  the  mass  rate 
of  water  accreted  (eq  6),  they  are  also  time-dependent 
in  the  same  manner  as  the  mass  rate  (i.e.  in  the  col- 
lection efficiency  E).  In  general,  then,  the  correct 
form  of  the  equations  exhibits  significant  time  de- 
pendence. A scheme  to  update  the  time-dependent 
parameters  as  ice  accretion  proceeds  is  described  later. 

In  summary,  the  major  changes  in  the  characteristic 
manner  in  which  ice  accretes  on  a given  object  with 
time  are:  the  variation  in  the  flow  field  (or  collection 
efficiency)  that  affects  the  subsequent  mass  rate  of 
collection,  an  increase  in  surface  area  that  affects  the 
total  heat  transfer  coefficient,  and,  to  a lesser  degree, 
changes  in  the  heat  capacity  of  the  object  through  the 
addition  of  ice. 

5.  Time  dependencies  of  lesser  order 

Ice  density 

The  effects  described  here  are  less  significant  but 
can  be  quantified  to  some  degree  and  may  be  worth- 
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while  including  in  a detailed  description  of  the  icing 
process.  The  first  of  these  is  the  ice  density.  Macklin 
(1962)  found  that  in  relatively  low  velocity  flows  and 
for  surface  temperature  below  0°C  (Case  3 above),  the 
ice  density  could  be  obtained  as  follows: 


(Mg/m3).  (16) 


for  ice  formed  either  near  the  melting  point  or  at  higher 
velocities,  or  for  larger  droplets,  the  density  will  ap- 
proach the  theoretical  maximum  density  of  0.91 7 Mg/m3 
and  eq  16  will  no  longer  apply.  This  variation  may  in- 
duce a time  dependence  since  the  volume  of  the  ice 
accretion,  as  well  as  its  mass,  will  change  the  subsequent 
flow  field.  The  volume  can  be  determined  by  the  mass 
deposited  (eq  6)  and  the  density  as  given  in  eq  16. 


Run-back  icing 

For  Case  2 (T  = 0),  only  a fraction  of  the  accreted 
water  is  initially  frozen.  In  our  first  formulation  we 
assumed  that  the  unfrozen  water  is  carried  away  by 
the  flow  after  imparting  its  supercooling  to  the  surface 
and  affects  neither  the  thermodynamic  processes  nor 
the  flow  field.  Observations  during  flight  tests  and  on 
experiments  have  indicated,  however,  that  icing  has 
occurred  on  the  rear  portions  of  objects  as  this  water 
flowed  back  over  the  surface  and  refroze.  This  forma- 
tion is  known  as  “run-back  icing."  It  is  a difficult 
parameter  to  include  because  of  the  formation  condi- 
tions-how  the  water  is  maintained  on  the  surface  and 
whether  heat  transfer  is  rapid  enough  to  freeze  it  be- 
fore it  is  partially  or  totally  removed  as  liquid  at  the 
trailing  edge.  Cansdale  and  McNaughton  (1977)  and 
Stallabrass  and  Lozowski  (pers.  comm.)  are  currently 
including  a run-back  formulation  to  determine  angular 
dependence  of  initial  icing  rate  without  time-dependence. 
It  appears  an  extension  of  the  present  work  with  this 
angular  dependent  thermodynamics  could  be  a signifi- 
cant step  forward  in  our  understanding  of  the  accretion 
process. 

Spongy  ice  growth 

Spongy  ice  is  formed  by  a mesh  of  interlocking 
dendrites  that  trap  a fraction  of  the  unfrozen  water 
and  cause  it  to  freeze  by  subsequent  heat  removal, 
possibly  through  conduction  processes  to  the  side  and 
rear  surfaces  of  the  object  on  which  ice  is  accreting. 
Spongy  ice  is  not  expected  to  occur  frequently  under 
flight  conditions  because  the  flow  would  mechanically 
disrupt  the  fragile  structure  that  this  ice  has  when  it 
first  begins  to  form.  However,  at  low  velocity  (<  1 0 
m/s,  List  1977)  or  low  Reynolds  number,  there  are 
indications  that  spongy  ice  would  exist. 


accreted  water).  E(W  - Wc)u  represents  the  initially 
unfrozen  fraction  of  the  total  accreted  water,  and  the 
curves  shown  give  the  fraction  lost  for  various  tempera- 
tures and  radii  of  the  objects.  The  quantity  (1 -fraction 
lost)  of  initially  unfrozen  water  is  that  which  contributes 
to  the  total  ice  accretion  through  the  spongy  ice  phenom- 
enon. 

In  the  numerical  model,  the  spongy  ice  growth  could 
be  incorporated  by  defining  a Reynolds  number  function 
that  determines  the  changeover  from  the  case  where  the 
unfrozen  fraction  is  shed  to  the  case  where  part  of  it  is 
incorporated  according  to  the  relations  shown  in  Figure 
4.  Experiments  are  necessary  to  verify  such  a condition 
since  the  thermodynamics  are  dependent  upon  geometry, 
the  heat  capacity  of  the  material,  and  the  characteristics 
of  the  flow  field  in  order  to  freeze  additional  water  over 
and  above  that  determined  for  Case  2 (7"  = 0)  and  may 
be  difficult  to  quantify. 

The  shape  of  the  accretion  is  especially  critical  since 
the  cross-sectional  diameter  (2 a)  can  be  modified  by  the 
spongy  growth  regime,  as  in  the  lobe  structures  seen  on 
hailstones.  These  changes  in  shape  lead  to  substantial 
changes  in  the  subsequent  mass  rate  of  accreted  water 
through  the  collection  efficiency  effect,  the  total  sur- 
face area  available  for  heat  transfer,  and  other  first-order 
feedback  effects  that  are  controlled  by  the  cross  section. 


NUMERICAL  ICE  ACCRETION  MODEL 


1.  Major  subroutines 

In  a first  attempt  to  include  time  dependence  and 
its  subsequent  effects  on  later  ice  accretion,  a numerical 
model  of  the  ice  accretion  processes  described  has  been 
programmed  on  a digital  computer.  The  model  is 
programmed  in  the  language  Basic  on  the  Dartmouth 
Time  Sharing  System  (DTSS),  which  uses  a Honeywell 
66/40  computer.  Hook-up  is  through  an  acoustic 
coupler  connected  by  telephone  line  to  the  computer. 
The  model  includes  software  using  DTSS  subroutines 
for  graphic  display  of  the  results  using  a Tektronix 
Model  4013  CRT  display.  A listing  of  the  program, 
called  Ice9,  is  given  in  Appendix  A.  In  Figure  5,  the 
feedback  effects  of  time  dependence  are  indicated  by 
the  various  arrows.  The  major  one  is  the  recomputation 
of  the  object’s  profile  after  ice  accretes,  based  on  the 
mass  of  accreted  ice  and  the  density  variation.  The  new 
profile  is  then  used  to  update  the  flow  field  and  deter- 
mines changes,  for  example  in  the  collection  efficiency 
E and  heat  transfer  coefficient  H,  for  the  next  time  step. 

Each  block  of  Figure  5 corresponds  to  a subroutine 
of  the  program,  and  we  will  describe  each  of  these, 
keeping  in  mind  their  overall  interconnection  as  shown 
in  Figure  S.  The  initial  variables  chosen  are  indicated 
in  Figure  6.  These  have  been  described  in  earlier 


Figure  4.  The  fraction  of  accreted  water  lost 
plotted  against  E(W  - WJ  u for  various  temper- 
atures and  average  radii  (cm)  of  accreting  ob- 
jects (after  Cams  and  Macklln  1973). 


In  a previous  paper  (Ackley  1977)  we  described 
how  the  spongy  growth  can  be  incorporated  into  a full 
time-dependent  solution  for  icing  prediction.  Here  we 
briefly  describe  the  variations  from  the  thermodynamics 
detailed  in  this  report. 

The  standard  thermodynamics  are  first  evaluated  to 
determine  that  the  T = 0°C  case  exists.  Instead  of  a 
frozen  fraction  F depending  only  on  the  available  heat 
transfer,  empirical  relations  are  used  to  determine  the 
fraction  of  the  accreted  water  that  is  lost  to  the  flow 
as  shown  in  Figure  4.  The  quantity  (W-WJu  is 
plotted  against  the  fraction  lost,  W and  representing 

the  total  twc  and  frozen  Iwc  respectively  (i.e.,  Wc  cor- 
responds to  that  amount  of  Iwc  that  would  raise  the 
surface  temperature  to  0°C  but  still  freeze  all  the 
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Figure  5.  Block  diagram  of  the  numerical  model  for  calculating  ice  accretion.  Each  block  represents  a subroutine  of 
the  program.  Flow  diagrams  and  pertinent  equations  for  each  of  these  are  given  in  Figures  6-9. 


Initialize  Variables 

Air  temperature  = T0 
Liquid  water  content  = Iwc 

Droplet  radius  = R0 

[or  droplet  radius  distribution  = )J 

Object  airspeed  = u 
[or  free  stream  air  velocity] 

Figure  6.  Initial  variables  input  in  the  run 
control  sequence  of  the  numerical  model. 
Other  parameters  such  as  the  object  size 
and  total  time  of  the  simulation  can  also 
be  varied  but  these  require  program 
statement  modifications  (see  Appendix  A). 


sections,  and  we  need  only  mention  here  that  the  pro- 
gram is  initially  set  for  a circular  cylinder  of  fixed 
radius  and  heat  capacity,  so  these  parameters  are  built 
in  as  program  statements  rather  than  initial  variables. 
They  can  be  varied,  however,  and  comments  within  the 
program  statements  (see  Appendix  A)  indicate  the 
appropriate  variables  that  refer  to  them.  The  choice 
of  droplet  sizes  is  controlled  by  a separate  subroutine, 
Dropdst,  which  will  be  described  later.  The  program  is 
initially  set  to  run  by  asking  a series  of  questions  con- 
cerning the  selection  of  initial  variables,  so  the  choice 
of  either  a single  droplet  size  or  a number  of  sizes 
(e.g.  mean  and  standard  deviation  for  a Gaussian  dis- 
tribution) is  automatically  routed  in  the  program. 

In  Figure  7a,  the  object’s  profile  is  recalculated  at 
a time  increment  set  for  the  given  icing  conditions. 

The  first  step  is  to  take  the  mass  of  the  icing  accretion 
output  by  the  thermodynamics  and  the  subroutine  for 
ice  density  (Fig.  7b)  and  compute  the  volume  of  the 


accreted  ice.  The  ice  density  depends  on  the  rela- 
tionship between  droplet  size,  impact  velocity  and  sur- 
face temperature  as  given  in  eq  16  (Macklin  1 962). 
Generally,  for  high  velocity,  high  surface  temperature 
cases  (>  -5°C)  the  ice  density  approaches  the  maxi- 
mum value  for  bubble-free  ice  of  0.91 7 Mg/m3.  The 
program  uses  eq  16  only  until  pice  reaches  0.917 
Mg/m3 ; the  density  then  remains  fixed  at  that  value.  The 
new  object  profile  is  chosen  to  conform  to  an  elliptical 
shape.  This  analytic  form  for  the  object  profile  con- 
siderably simplifies  the  recomputation  of  the  flow 
field,  and  more  elegant  solutions  from  potential  theory 
can  be  applied.  Whether  the  forms  chosen  are  in  fact 
correct  is  the  subject  of  experimental  investigation 
(Acklty  et  al.  1978)  that  will  provide  the  basis  for  a 
better  estimate  of  the  object  profile.  At  present,  the 
minor  axis  of  the  ellipse  perpendicular  to  the  flow 
direction  remains  constant  at  the  initial  value  while  the 
major  axis  parallel  to  the  flow  direction  is  increased  in 
the  direction  of  the  accreting  ice  in  relation  to  the 
mass  and  density  of  the  accreted  ice  (e.g.  see  Fig.  1 1, 
13, 14).  Preliminary  experiments  indicate  that  these 
assumptions  are  not  seriously  in  error,  although  ob- 
served small  changes  in  cross  section  can  markedly 
affect  the  results,  since  the  cross  section  appears 
directly  in  both  the  mass  rate  and  the  thermodynamic 
equations. 

Figure  8 shows  the  subroutines  controlling  the 
computation  of  the  droplet  trajectory,  the  collection 
efficiency  and  trajectory  control,  and  the  free  stream 
velocity.  The  subroutine  Traject  (Fig.  8a)  integrates 
the  equation  of  motion  for  the  droplet,  calling  upon 
the  subroutine  Velocity  (Fig.  8c)  as  needed  to  deter- 
mine the  velocity  of  the  flow  field  at  each  point  on 
the  'roplet’s  path,  and  updates  the  droplet’s  position 
ur  ,,t  one  of  the  three  cases  of  interaction  with  the  ob- 
ject is  determined.  In  order  to  compute  the  collection 
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Droplet  trajectory  calculation 

[TRAJECT] 


J 


f 


1 Calculate  k,  (inertial  parameter)  k = \/a  = 2pw  /?drop  ul9r>a 
(Langmuir  and  Blodgett,  1946) 

2 Calculate  Reynolds  number  of  the  droplet  moving  with  respect  to  the  air  -4 

3 Calculate  the  drag  coefficient  cD(cD  Re/24  = fl6mtudRdrop): 

Beard  and  Pruppacher  (1969)  give 

cD  = 1 + 0.102  Re-95*  0.2  < Re  < 2 

cD  = 1 + 0.1 1 5 Re-802  2 < Re  < 21 

cD  = 1 +0.189  Re-632  21  < Re  < 200 

4 Call:  subprogram  VELOCITY  for  calculating  free  stream  velocity  of  air 
at  droplet’s  current  location 

5 a)  dud  = cD  Re/24  (ud  - u0)  (dt/k)  (change  in  velocity  of  droplet) 

b)  ud  = ud  + dud  (new  velocity  of  droplet) 

c)  xd  = xd  + ud  dt  (new  position  of  droplet) 

6 Is  the  droplet  downstream  far  enough  for  it  to  be  capable  of  hitting  the 
profile,  i.e.  is 


1 


Xrf  < a 


No 


Reiterate 


7  Compare  the  tangent  profile  and  the  droplet  trajectory  tangent  (v0/u0) 

if  I ^profile  I < I (*V"o)  I No 

yes,  then  particle  (could  have)  missed  profile 


8  Check  to  see  by  how  much  droplet  missed  profile 

if  /^+jd)_1<.oi 

\ a2  b2  / 


then  droplet  missed  profile,  but  the  trajectory  was  •«  tangent  

9 **  missed  profile  **  let /V0  = 1 

end  trajectory 

10  Did  we  hit  target  yet?  /— V + - 1 < 0 

W \bf  No  Reiterate  — 

1 1 Was  the  trajectory  within  .1  arc  tan  units  of  being  tangent  to  the  profile 

No  

12  **  tangency  found  **  Let  HO  = 2 

end  trajectory 

13  Let  HO  = 3 -4 

*♦  incorrect  hit  ** 


a.  Subroutine  Tra/ect  uses  the  object  profile  output  (Fig.  7a)  and  calls  the  subroutine 
Velocity  to  calculate  the  free  stream  velocity  (Fig.  8c).  It  then  Integrates  the  equation 
of  motion  to  determine  whether  and  how  a droplet  of  particular  size  Impacts  the  ob- 
ject. The  HO  values  are  used  to  Indicate  the  droplet  has  missed  the  object  (HO  = 1), 
hit  the  object  at  a tangent  (HO-  2),  or  hit  the  object  but  not  at  a tangent  (HO*  3). 

Figure  8.  Subroutines  controlling  computation  of  the  droplet  trajectory,  collection 
efficiency  and  trajectory  control,  and  free  stream  velocity. 
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Collection  efficiency  calculation  and  trajectory  control 

[COLL  EFF] 


TRAJECT 


b.  Subroutine  Coll  Eff  takes  the  droplet  trajectory  Information  to  compute  the  collection  efficiency  E for  the 
given  droplet  size  and  object  profile.  If  the  correct  tangency  value  of  droplet  impact  is  not  found,  the  sub- 
routine readjusts  the  droplet's  Initial  off-axis  position  (y$)  and  calls  Tra/ect  to  compute  the  correct  tangent 
trajectory.  Once  the  collection  efficiency  is  assigned,  the  program  proceeds  to  Thermo  (Fig.  9a). 


Figure  8 (cont'd). 


) 

I 
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temperature.  Then,  based  on  the  current  surface  tem- 
perature, the  appropriate  thermodynamic  expression 
is  chosen  (Fig.  9a).  The  complete  thermodynamic 
balance  yields  the  new  surface  temperature,  or  new 
ice  fraction  accreted  if  in  the  T = 0°C  temperature 
regime.  Of  course,  care  must  be  taken  to  maintain 
physical  as  well  as  thermodynamic  continuity  as  the 
temperature  proceeds  through  the  limits  of  each 
thermodynamic  argument.  The  program  checks  for 
errors  in  continuity  and  redirects  execution  at  the 
appropriate  level. 

At  the  end  of  the  icing  thermodynamics,  the  output 
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is  given  as  new  ice  mass,  total  heat  capacity,  and  surface 
temperature.  After  several  iterations  of  the  thermo- 
dynamics, the  new  total  ice  mass  is  fed  back  to  the  ob- 
ject profile  calculation  (Fig.  5,  7)  to  recompute  the 
object’s  size,  flow  field  characteristics,  and  surface  area 
for  the  next  time  period  of  the  calculation. 

2.  Options  for  droplet  size  variations 

In  addition  to  the  basic  computational  portions  just 
discussed,  the  Dropdst  subroutine  allows  the  specifica- 
tion of  more  than  one  droplet  size  for  any  given 
simulation.  Dropdst  has  three  options  available  on 


Stream  velocity 

(VELOCITY] 


From  complex  potential  flow  theory  the  velocity  at  any  point  x,  y in  the 
flow  around  an  elliptical  cylinder  is  given  by: 


where: 

a=  ( (x2  -y1  - c2)2  + 4x2y2\  ^ 


R = ( x 2 +y2)'/2  = xd 


0 = tan*1  (y/x ) 
a = semimajor  axis  length  (||  u) 
b = semiminor  axis  length  (1«) 


see  Mllne-Thomson,  1 960 

c.  Subroutine  Velocity  computes  the  velocity  In  the  flow  at  a given  x,y 
position  based  on  potential  flow  theory.  Traject  (Fig.  8a)  uses  this  In- 
formation in  Integrating  the  equation  of  motion  and  computing  the 
droplet's  trajectory. 

Figure  8 (cont'd).  Subroutines  controlling  computation  of  the  droplet 
trajectory,  collection  efficiency  and  trajectory  control,  and  free  stream 
velocity. 


droplet  distributions  in  addition  to  the  single  droplet 
size,  which  does  not  require  the  calling  of  the  Dropdst 
subroutine.  Currently  the  droplet  distributions  avail- 
able are  a Gaussian,  an  Erlang  of  order  1 , 2 or  3,  or  a 
histogram  of  sizes  that  could  be,  for  example,  from 
experimental  data.  In  running  the  program,  the 
operator  is  asked  to  specify  whether  a distribution 
is  desired  and,  if  so,  this  subroutine  is  called.  It  then 
asks  additional  questions  concerning  the  type  of  dis- 
tribution desired  and  the  necessary  parameters  for 
specification  of  this  distribution. 

Gaussian  distribution 

The  program  asks  for  the  mean,  radius  variance  and 
total  number  of  droplet  sizes  desired  up  to  a maximum 


of  nine  sizes.  The  subroutine  then  apportions  the 
available  liquid  water  content  (specified  initially  in  the 
control  sequence)  into  these  categories  in  such  a way 
that  the  total  number  of  droplets  conforms  to  the 
distribution  as  specified.  The  Gaussian  is  sampled 
between  ± 2 standard  deviations  of  the  mean  so  that 
it  includes  95%  of  the  distribution. 

Erlang  distribution 

The  Erlang  distribution  is  of  the  form 


Fr(r-r0) 


X*  ( r - r0)*_1  e-*(r-r<)) 

(srnri 


(17) 
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Icing  thermodynamics 

[THERMO] 


1 Calculate  the  ambient  water  vapor  pressure  (£0) 

Smithsonian  Meteorological  Tables,  1951 

2 Call:  subprogram  HTTRANS1  for  calculating  heat  transfer  coefficient 

(watts/°C  (meter  of  blade/2)  sec| 

3 Calculate  mass  of  water  collected/meter  of  blade  - sec 

m = £ • Iwc  ■ u ■ b 


4 r < o 

branch  to  appropriate  thermodynamics 

initially,  , 

T = 0 

(see  E.A.  Brun,  1957) 

T _ r°  + 3o<5ff 

T > 0 

i.e.  recovery  temp 

S T > 0 what  is  0? 
a)  Q = ~h\t-T[ 


» [r  ■ - To  - ^ ♦ .6  L0(E'  - £0)/p  cp  ^ - 


R - surface  recovery  factor 
L0  - latent  heat  of  vaporization  of  water  @ 0°C 
£'  = vapor  pressure  of  water  at  surface  temperature 

b)  Calculate  new  T 

T ~ (^cpmblade  + ^)/£pmbladc 

c)  Is  T still  above  0°C? 

• * end  THERMO  •* 

6 T = 0 


+ 


.6  Lq(E  - Eq)Ip  Cp  aj^J  - 

-mcpHtO  ( - T0 

\ cp  HjO  “pH, O' 


£ * fraction  of  collected  water  accreted  as  ice 
L - latent  heat  of  freezing 

b)  £ = F-QILm 

c)  is  £ > 0 Temperature  above  0°C  No  ■ 

d)  is  £ < 1 Temperature  goes  below  0°C  No  ■ 

*)  cpm  blade  = cpm  blade  +m^cplce 

*•  end  THERMO  •* 

7 T < 0 


T0  - * •6tot£  “fo)/£cpalrJ  “ 


a.  Subroutine  Thermo  calculates 
mass  of  water  arriving  and  con- 
verted Into  Ice  or  shed  as  un- 
frozen water  if  T > 0°C. 

Figure  9.  Based  on  the  ambient 
conditions,  the  calculation  of  mass 
accreted  is  completed  by  Incorp- 
orating the  thermodynamics  at 
the  surface. 


[rieisi-  - r( 
V CPH,0 


~L 

cpH,0  “pH,0/ 


b)  T * (T’fpmblide  + Wcpmblade 
e)  cpm blade  = cpm blade  + mcple* 

SUBEND 
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Heat  transfer  coefficient 

[HTTRANS1 1 


1 H0  = 31.01  u°0S6l(2a)M  /^imiiV56  .watts  . 

\ 1014  F*  j m^Csec 

F*  = film  temperature  = ^ 

2 the  total  heat  transfer  coefficient/111—^  ^lade 

b.  Subroutine  Httransl  computes  the  total  heat  transfer  coef- 
ficient based  on  the  current  surface  temperature  and  object 
dimensions.  It  Is  called  by  Thermo  (Fig.  9a)  as  needed. 

Figure  9 (cont ’d).  Based  on  the  ambient  conditions,  the  cal- 
culation of  mass  accreted  Is  completed  by  Incorporating  the 
thermodynamics  at  the  surface. 


Figure  1 0.  Forms  of  the  Erlang  distribution  for  various  values  of  the 
parameter  k as  shown.  The  "start  radius  "r0  Is  zero  for  these  cases. 


and  has  the  forms  shown  in  Figure  10  for  various 
values  of  the  "Erlang  order” k and  for  “start”  values 
r0  = 0.  The  mean  or  expected  value  of  the  distribution 
is 

T=$+r0.  (18) 

As  Figure  1 0 shows,  a variety  of  distributions  can 
be  obtained  by  specification  of  a few  parameters 


(h,  X,  r0),  making  it  easy  to  change  the  characteristics 
of  the  distribution  for  simulation  comparisons.  The 
distributions  are  also  more  reminiscent  of  the  variety 
one  sees  in  real  cloud  droplet  distributions  than  the 
Gaussian  distribution  (Mason  1971).  The  subroutine 
asks  the  specification  of  the  values  of  k,  T,  r0  and  the 
number  of  droplet  classes,  up  to  a maximum  of  nine, 
that  are  desired.  The  subroutine  then  calculates  the 
value  of  X from  eq  18  and  apportions  the  available  Iwc 
into  the  specified  classes  according  to  eq  1 7.  The 
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Figure  1 1.  (Top)  The  change  in  profile  dimension  at  50 -s  intervals  by  ice 
accretion  is  indicated  by  the  profile  shapes  on  the  Initial  half  cylinder.  The 
black  line  coming  from  the  right  Is  the  tangent  trajectory  at  the  beginning 
(top  of  the  black  line)  and, end  (bottom  of  the  black  line)  of  the  Icing  period 
of 300  s.  (Bottom)  The  surface  temperature  of  the  front  half  cylinder  as  a 
function  of  time  Is  plotted.  As  seen,  the  surface  temperature  reaches  an 
equilibrium  value  within  about  30  s after  the  Icing  starts. 


subroutine  samples  the  distribution  from  the  start  radius 
r0  to  the  value  r0+2xr 

Histogram  distribution 

Dropdst  also  allows  the  specification  of  an  arbitrary 
distribution  such  as  may  be  obtained  during  experi- 
ments or  in  airborne  sampling.  In  this  case,  the  input 
to  the  subroutine  is  the  number  of  droplet  sizes  (up  to 
nine)  and  pairs  of  values  that  indicate  the  droplet  radius 
and  its  percentage  contributions  to  the  total  numbers. 
The  subroutine  then  computes  the  volume  percentages 
and  the  total  number  required  to  equal  the  designated 
Iwc. 

3.  Option  for  the  helicopter  rotor  case 

As  mentioned  in  the  introduction,  helicopters  have 
serious  icing  problems  since  they  normally  operate  at 
the  altitudes  where  icing  conditions  most  often  occur. 
They  also  have  the  interesting  characteristic  that  the 
linear  velocity  of  the  rotating  blade  varies  over  a large 
range  (eq  1 ),  from  essentially  near  zero  in  the  hub 
region  to  nearly  sonic  at  the  blade  tip.  Interest  in  how 
these  changing  conditions  influence  ice  accretion  pro- 
cesses was  a continuing  impetus  in  the  development  of 
this  numerical  scheme.  Two  options  to  include  the 
velocity  variability  were  therefore  built  into  the  pro- 
gram. 

The  first  of  these  selects  specific  values  of  the 
velocities  and  computes  the  ice  mass  and  profile  changes 


as  a function  of  time.  Figure  1 1 shows  such  a calcula- 
tion for  a velocity  of  10  m/s  and  initial  conditions  of 
0.00143  kg/m3  Iwc,  droplet  radius  of  45  pm  and 
Ta  = -1 7.7°C.  Any  number  of  these  velocities  can 
be  selected,  allowing  details  of  the  ice  accretion  process 
to  be  seen,  such  as  the  time  dependence  of  the  collection 
efficiency  or  heat  transfer.  The  heavy  black  line  is  the 
tangent  trajectory  used  to  compute  the  collection 
efficiency  for  this  case.  Examples  of  these  will  be 
shown  in  the  next  section.  The  figure  illustrates  the 
three  main  time  intervals  of  the  calculation.  The 
spacing  of  the  dots  in  the  temperature  plot  indicates 
the  time  period  (2  s)  that  we  take  for  the  pseudo  steady- 
state  heat  balance,  i.e.  in  the  2-s  time  interval  the 
thermodynamic  properties  of  the  system  do  not  change. 
After  25  repetitions  of  the  thermodynamic  balance  with 
successive  updates  of  the  heat  flux,  the  profile  dimensions 
are  sufficiently  changed  to  require  the  updating  of  the 
collection  efficiency,  ice  mass  and  profile  shape  charac- 
teristics. These  updates  take  place  at  50-s  intervals,  as 
indicated  by  the  hacks  on  the  time  scale.  The  changes 
in  the  profile  dimensions  indicated  on  the  figure  take 
place  at  these  intervals.  The  third  interval  is  the  total 
time  of  the  calculation,  indicated  by  the  full-scale 
length  of  the  time  scale  (300  s).  Programming  mechanics 
dictate  that  these  times  be  integer  multiples  of  the 
smallest  value.  However,  they  can  be  modified  within 
this  constraint. 
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Model  time:  300  sec 


Ta  = - 15.5  (°C) 

Iwc  = 0.00108  kg/m3 


Figure  12.  Numerical  simulation  of  helicopter  rotor  blade  Icing 
showing  the  surface  temperature,  leading  edge  ice  thickness, 
collection  efficiency  and  fraction  accreted  (as  Ice)  as  a function 
of  velocity  for  the  Initial  conditions  as  shown.  The  maximum 
Ice  thickness  occurs  when  heat  transfer  conditions  optimize 
with  amount  of  accreted  water  (eq  6)  to  form  the  most  ice. 


The  second  display  option  is  shown  in  Figure  12. 
Here  the  three  quantities  equilibrium  surface  tempera- 
ture, leading  edge  (maximum)  ice  thickness,  and 
variations  of  collection  efficiency  and  fraction  of  col- 
lected water  accreted  as  ice  are  given  as  a function  of 
velocity  along  the  blade  from  0 m/s  to  250  m/s.  The 
program  computes  and  stores  the  heat  balance,  ice 
thickness,  and  collection  efficiency  for  each  10-m/s 
velocity  interval  for  a period  of  300  s and  then  steps 
to  the  next  velocity  value.  When  this  has  been  accom- 
plished for  25  velocity  values  (0  to  250  m/s  in  10-m/s 
intervals),  the  information  saved  at  the  end  of  each 
calculational  step  is  displayed.  Straight  lines  are  drawn 
between  the  points  to  complete  the  curves  as  shown. 

In  addition,  the  lowest  plot  shows  the  percentage  of 
the  accreted  water  that  is  frozen  into  ice.  As  shown 
for  this  particular  plot,  the  fraction  accreted  and  frozen 
drops  below  100%  at  the  point  where  the  surface  tem- 
perature rises  to  0°C  as  the  thermodynamic  equations 
would  indicate. 

As  shown  in  Figures  1 1 and  1 2 the  program  also 
includes  software  for  plotting  data.  These  are  given 
by  the  subroutines  Plot  T,  Pt  Blade,  Plot  D,  Plot  M 
and  Thermopt.  The  function  of  each  of  these  is 
described  in  the  comments  accompanying  the  program 


listing  in  Appendix  A.  In  addition  to  the  plots  shown 
in  Figures  1 1 and  1 2,  the  program  can  also  show  a plot 
of  the  droplet  size  distribution  if  a distribution  is 
selected.  The  variation  of  the  thermodynamic  quantities 
given  in  eq  7,  10,  1 1 and  12  as  a function  of  segment 
velocity  (similar  to  the  surface  temperature  plot  in 
Figure  12)  can  also  be  plotted.  The  quantities  plotted 
are  the  convective  flux,  evaporative  flux,  kinetic  energy 
of  impact  term,  the  latent  heat  requirement  to  freeze 
the  accreted  water,  and  the  influx  term  of  heat  (nega- 
tive) due  to  droplet  supercooling.  These  are  defined  in 
eq  7,  10, 11  and  12  and  discussed  more  fully  in  the 
following  section. 

RESULTS 

In  this  section,  we  give  some  examples  of  the  effect 
of  including  the  time  dependence  as  formulated  pre- 
viously. First,  the  effects  of  including  collection  ef- 
ficiency feedback  and  a distribution  of  droplet  sizes 
are  shown  in  Figure  13. 

Here  the  amount  of  ice  accreted  as  time  proceeds  is 
plotted  for  the  initial  conditions  of  Iwc  = 0.00108 
kg/m3,  initial  ambient  air  temperature  = -5.5°C,  and 
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Figure  13.  Droplet  trajectories  and  ice  profile  changes  for  a Gaussian  droplet 
size  distribution.  Droplet  sizes  and  respective  liquid  water  contents  are  given 
in  the  inset;  the  collection  efficiencies  for  the  various  droplet  categories  and 
their  changes  with  profile  dimension  changes  are  given  on  the  right  hand  scale. 
Profile  Ice  thickness  changes  are  shown  at  50-s  intervals.  Thermodynamic 
conditions  for  this  simulation  are  such  that  the  surface  temperature  was  at 
0°C  (equilibrium)  and  some  accreted  water  was  not  frozen. 


velocity  = 60  m/s.  The  droplet  distribution  selected 
is  a Gaussian  with  six  classes.  The  mean  radius  is  1 1pm 
(17x  10-6  m)  and  the  droplet  radius  standard  deviation 
is  8pm  about  the  mean.  As  shown  by  the  inset  table 
listing  the  Iwc  amounts  by  droplet  size,  very  little  of 
the  total  water  available  is  partitioned  into  the  lower 
droplet  sizes  (for  the  two  lowest  droplet  sizes  the  Iwc 
is  2 to  3 orders  of  magnitude  lower  than  for  any  of  the 
four  higher  classes). 

Referring  to  the  profile  plots,  the  dark  lines  origi- 
nating at  the  right  refer  to  the  collection  efficiency 
and  its  changes  with  time  and  droplet  size  as  shown. 
The  thickness  of  these  lines  results  from  "overprinting” 
on  the  display  screen  so  the  top  of  the  line  for  any 
droplet  class  is  the  tangent  trajectory  that,  when 
scaled  to  the  object  radius,  gives  the  initial  value  of 
collection  efficiency  for  that  class  (scale  at  right).  The 
bottom  of  each  thick  line  is  the  final  collection  ef- 
ficiency value  after  300  s of  ice  accretion  has  taken 
place.  For  the  largest  droplets  (>  20.2pm)  the  col- 
lection efficiencies  are  generally  0.75  or  above  for 
this  velocity.  The  collection  efficiency  drops  markedly 
for  sizes  below  20pm  and  the  decrease  in  collection 
efficiency  with  ice  accretion  (given  by  the  distance 
between  the  top  and  bottom  trajectories  or  thickness 
of  each  line)  is  also  greater  for  the  smaller  droplets. 

The  ice  accretion  rate,  shown  by  the  changes  in  the 
object’s  profile  at  50-s  intervals,  does  not  change 
dramatically  since  90%  of  the  liquid  water  is  con- 
centrated in  the  three  largest  droplet  'lasses.  The 
collection  efficiencies  for  these  classes  v^ry  by  the 
smallest  amounts  (~  3 to  5%)  at  this  velocity  and 


object  radius,  so  the  rate  of  ice  thickness  change  is  not 
substantial  for  the  time  period  of  this  simulation.  It  is 
easily  seen  from  this  plot,  however,  that  a lowering  of 
the  mean  droplet  size  such  that  a significant  shift  occurs 
in  the  amount  of  Iwc  available  in  small  droplets  (<14pm 
radius)  would  substantially  change  the  ice  accretion  rate. 
The  collection  efficiency  and  its  change  with  object 
profile  are  both  more  markedly  affected  at  the  lower 
droplet  sizes;  for  example  the  collection  efficiency  of 
the  1 -pm  category  changes  from  1 5%  to  6%,  a 60% 
change  with  the  object  profile  change. 

We  illustrate  this  further  in  Figure  14  where  two 
profiles  are  shown  for  the  same  conditions  except  that 
the  droplet  size  is  given  as  30pm  on  the  top,  and  half 
that,  15pm,  on  the  bottom.  From  this  figure  we  see 
that  the  total  center  line  ice  thickness  is  reduced  by  40% 
by  changing  the  droplet  radius,  even  though  the  ambient 
temperature  and  liquid  water  content  are  the  same  for 
both  cases.  Similarly,  the  calculated  surface  tempera- 
ture is  higher  for  the  larger  droplets  since  more  latent 
heat  is  liberated  per  unit  time,  which  increases  the  heat 
flux  to  the  object  surface  and  increases  its  temperature. 

A more  complete  treatment  of  the  comparison  be- 
ween  experimental  results  and  model  simulations  will 
be  the  subject  of  later  reports.  These  comparisons  are 
important  because  of  their  use  in  updating  model  re- 
sults and  improving  on  portions  such  as  the  object  pro- 
file shape  and  heat  transfer  coefficient  which  are 
currently  based  on  assumptions  (e.g.  the  time  period 
of  steady-state  thermodynamics)  that  may  require 
modification.  In  Figure  1 5,  we  compare  and  discuss 
one  such  simulation  with  experimental  data  to 
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Iwc  - 0.00143  kg/m3 
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Figure  14.  Accreted  ice  thickness  and  temperature  for  the  ambient  con- 
ditions indicated.  The  profile  updates  occur  at  50*  intervals.  The  differ- 
ence in  accreted  ice  thickness  for  the  two  cases  occurs  because  of  the 
change  in  collection  efficiency  for  the  droplet  sizes  30pm  (top)  and  15pm 
(bottom). 


indicate  both  the  general  validity  of  the  model  and 
some  specific  differences  that  require  more  experi- 
mentation at.d  better  model  parameterization.  This 
figure,  in  the  ice  thickness  distribution  portion,  shows 
the  simulated  and  experimental  results.  This  experi- 
mental result  has  been  previously  discussed  in  Ackley 
(1977)  and  is  more  fully  described  there.  Briefly,  the 
measurements  were  taken  from  a rotating  cylinder 


system  at  3600  rpm  in  a coldroom  at  the  ambient 
temperature  shown.  Water  of  the  mean  droplet  size 
and  Iwc  indicated  was  sprayed  onto  the  rotor.  Com- 
paring the  two  ice  thickness  curves  it  is  seen  that  the 
simulated  ice  thickness  is  somewhat  less  than  the 
experimental  accumulation  but  is  reasonable  within 
the  constraints  imposed  by  both  the  assumptions  in 
the  simulation  (e.g.  one  droplet  size  at  the  mean 
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Figure  15.  Comparison  between  experimental  data  (ambient  con- 
ditions inset)  and  model  simulation  of  the  same  conditions.  Agree- 
ment is  reasonable  considering  experimental  errors  in  Iwc  (-  25%) 
and  use  of  a single  droplet  size  (mean  experimental  size)  in  the 
simulation  compared  to  a distribution  in  the  experiment. 


experimental  value)  and  the  errors  in  measurement  in 
the  experiment  on  both  droplet  size  and  Iwc  (-25%). 

In  the  experiment,  analysis  of  the  ice  properties 
showed  a grain  size  change,  indicating  a transition  in 
surface  temperature  from  T < 0°C  to  T = 0°C  in  the 
vicinity  of  the  blade  where  the  line  . velocity  was 
20  m/s.  As  shown  by  the  calculated  temperature 
distribution  (top  of  Fig.  15)  the  calculated  rise  to  0°C 
takes  place  in  about  the  same  velocity  region.  This 
agreement  would  indicate  the  simulated  heat  transfer 
and  thermodynamic  relations  are  approximately 
compatible  with  the  experimental  measurements. 

The  location  of  the  maximum  ice  thickness  in  the 
experimental  data  differs  from  that  in  the  simulation. 

The  latter  indicates  that  the  ice  thickness  should  be 
increasing  all  the  way  to  the  tip  of  the  experimental 
rotor  (v~  100  m/s)  for  this  set  of  conditions.  How- 
ever, a set  of  experiments  in  which  the  Iwc  and  droplet 
distributions  were  held  constant  while  the  ambient 
temperature  was  varied  showed  that  the  location  of 
the  experimental  maximum  ice  thickness  varied  only 
slightly  with  ambient  temperature  at  roughly  the  same 
location  shown  in  Figure  1 5.  This  continuity  be- 
tween the  experiments  indicates  the  maximum  position 
may  be  controlled  by,  for  example,  a flow  field  perturba- 
tion induced  by  the  blunt  end  of  the  cylindrical  rotor. 

We  plan  further  experiments  at  different  rotation  rates 
(i.e.  changing  the  velocity  distribution  without  changing 
the  physical  dimensions  of  the  rotor)  to  see  if  and  how 


this  effect  is  modified.  At  this  time,  therefore,  the 
difference  in  the  position  of  the  thickness  maximum 
between  the  experiments  and  the  simulation  cannot 
be  resolved  until  further  experimental  results  are 
available  that  can  test  the  end  effects  of  the  rotor 
tip. 

In  Figure  12  we  showed  the  accumulated  ice  for 
the  conditions  as  indicated.  Figure  16  gives  the 
relative  contributions  of  the  thermodynamic  terms 
for  the  same  conditions,  which  we  now  discuss  more 
fully.  The  heat  balance  requires  that  the  algebraic 
sum  of  the  positive  and  negative  contributions  to  the 
heat  flux  equal  zero  as  long  as  ice  is  accreting  (Fig.  16). 
If  the  balance  goes  positive,  then  the  excess  heat  is 
accounted  for  by  a rise  in  surface  temperature  (eq  14). 
At  velocities  below  50  m/s,  Figure  16  shows  that  the 
convective  flux  dominates,  contributing  about  50%  of 
the  necessary  heat  flux  to  freeze  the  accreted  water, 
with  smaller  and  nearly  equal  contributions  from 
droplet  supercooling  and  evaporative  heat  flux  of 
approximately  25%  each.  At  about  80  m/s,  the  con- 
vective flux  is  at  its  maximum  negative  value;  it  then 
becomes  less  negative  with  increasing  velocity  because 
of  the  contribution  to  this  flux  by  adiabatic  com- 
pressional  heating  of  the  flow  (/?w$/2cp  ajr  term  in 
eq  7).  Above  80  m/s,  an  additional  positive  term,  the 
kinetic  energy  imparted  by  droplet  impact,  also  be- 
comes significant.  The  increasing  positive  contribution 
of  these  terms  causes  the  amount  of  heat  to  balance 
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Figure  16.  Magnitude  of  the  individual  thermodynamic  terms  as  a 
function  of  velocity  for  the  simulation  shown  in  Figure  12.  The  terms 
are  defined  in  Figure  9a  and  by  eq  7,  10,  11  and  12. 


the  latent  heat  of  freezing  to  peak  and  then  drop  off 
since  the  other  negative  terms,  the  evaporative  and 
supercooling  influx,  decrease  more  slowly  with 
velocity  than  the  positive  terms  are  increasing.  The 
maximum  available  latent  heat  requirement  corre- 
sponds to  the  maximum  ice  thickness  at  140  m/s 
since  we  are  in  the  0°C  surface  temperature  region 
(Fig.  12,  top).  At  the  maximum  ice  thickness,  the 
heat  necessary  to  freeze  the  accreted  ice  is  about  45% 
from  droplet  supercooling,  25%  from  evaporative  flux 
and  30%  from  convective  flux.  Kinetic  energy  of  the 
droplets  consumes  about  7%  of  the  heat  at  this  point, 
the  remainder  being  balanced  by  the  latent  heat  of 
freezing.  No  heat  is  available  for  freezing  from  con- 
vective flux  above  200  m/s  due  to  the  heating  effects 
of  adiabatic  compression.  At  about  240  m/s,  the 
kinetic  energy  of  droplet  impact,  the  now  positive 
convective  flux,  and  the  latent  heat  of  freezing  all 
consume  equal  amounts  of  the  negative  flux  available 
from  evaporative  flux  and  droplet  supercooling. 

CONCLUSIONS  AND  FUTURE  STUDIES 

As  pointed  out  earlier,  one  important  result  of  the 
numerical  model  is  that  a more  accurate  comparison 
of  theory  with  experimental  results  is  made  possible 
by  including  time  dependence  in  a measurable  way. 


These  comparisons  can  provide  more  confidence  in 
the  simulation  of  conditions  that  are  not  experimentally 
known.  For  example,  it  is  clear  from  Figures  13  and 
14  that  a change  in  either  mean  droplet  size  or  in  the 
distribution  such  that  large  numbers  of  small  droplets 
compose  much  of  the  liquid  water  content  can  radically 
change  the  amount  of  ice  accreted.  Using  ground  icing 
sprayers  with  large  droplet  sizes  for  tests  representing 
cloud  conditions  may  therefore  be  highly  inaccurate. 

As  mass  accretes  the  collection  efficiency  changes  with 
time  and  is  also  a function  of  droplet  size.  The  numer- 
ical solution  graphically  portrays  (Fig.  13,  14)  how 
these  conditions  can  affect  the  ice  buildup  for  a given 
set  of  initial  conditions. 

In  a situation  such  as  icing  of  a helicopter  rotor 
blade,  the  thermodynamics  determining  the  ice  buildup 
at  any  particular  linear  velocity  can  show  wide 
variability  as  to  which  terms  dominate,  so  a numerical 
solution  and  graphic  display  can  provide  a convenient 
method  of  seeing  this  variability  (Fig.  12,  16).  For 
design  purposes,  additional  terms,  for  example,  in- 
creasing the  heat  flux  to  the  blade  surface  by  electrical 
heating,  can  be  evaluated  to  determine  whether  particu- 
lar icing  conditions  would  be  relieved  or  unaffected  by 
the  extra  heat  flux.  An  optimization  procedure  could 
then  oe  used  with  statistics  on  meteorological  conditions 
to  evaluate  the  likelihood  of  the  design  change  signifi- 
cantly increasing  the  chances  for  successful  operation. 
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Future  changes  in  the  model  will  include  a continu- 
ing comparison  with  experimental  data  to  determine 
more  exactly  the  physics  controlling  the  icing  process. 
Evidence  exists  that  cross  section  and  object  shape  can 
change  considerably  with  changing  meteorological 
conditions,  especially  at  higher  temperature  (Dickey 
1952),  so  a change  in  the  way  these  parameters  are 
included  is  probably  necessary.  Preliminary  compari- 
sons with  experiments  at  temperatures  below  -25°C 
(Ackley  et  al.  1978)  indicate  that  the  observed  total 
ice  thickness  after  a few  minutes  of  accretion  is  less 
than  what  is  predicted  in  the  current  simulation.  This 
difference  may  result  from  a more  rapid  lowering  of 
collection  efficiency  with  ice  accretion  than  that  cur- 
rently used,  necessitating  a reevaluation  of  the  flow 
field  taking  into  account  that  the  object  shape  is  irregu- 
lar and  not  elliptical.  List  (1977)  has  also  indicated  that 
use  of  potential  flow  should  be  modified  to  include 
significant  boundary  layer  effects  near  the  object. 

These  are  not  now  in  ihe  model,  but  can  be  parameter- 
ized as  long  as  experimental  evidence  indicates  they 
are  necessary.  Finally,  the  most  difficult  parameter- 
ization involves  the  thermodynamics  since  it  is  the 
most  complex  and  relies  on  the  most  assumptions. 

A major  reliance  is  on  an  assumed  heat  transfer  coef- 
ficient based  on  our  best  present  knowledge.  This 
coefficient  is  then  used  to  derive  the  highly  important 
surface  temperature.  Experiments  using  isotopic 
analysis  (deuterium:  hydrogen  and  oxygen-1 8:oxygen- 
16  ratios)  on  the  accreted  ice  can  perhaps  be  used  to 
derive  the  surface  temperature  and  independently 
verify  our  choice  of  heat  transfer  coefficient.  The 
importance  of  heat  flow  through  the  back  surface  of 
the  object  and  the  variation  of  the  heat  transfer  with 
time  and  conditions  may  then  be  better  parameterized 
for  model  prediction  use. 

At  present,  within  the  data  limitations  that  have 
been  used  for  its  formulation,  the  model  offers  a 
significant  method  of  evaluating  time  dependence  (or 
independence)  of  the  various  parameters  that  influence 
the  icing  process.  With  this  information,  extension 
to  engineering  design  may  be  facilitated  by  easy  access 
to  a number  of  simulated  “case  histories"  applying  to 
aircraft  icing  or  ground-based  structural  problems. 

Total  icing  loads  can  be  computed  by  breaking  the 
structure  down  geometrically,  treating  each  shape 
individually,  and  summing  up  the  individual  loads 
into  the  total  for  the  structure. 

The  beauty  of  the  model  lies  in  its  subroutine 
construction,  which  allows  updating  of  particular 
aspects  (flow  field,  thermodynamics)  as  more  evidence 
becomes  available,  without  the  need  for  completely 
altering  those  portions  that  are  substantially  proven. 


The  utility  of  including  a more  complex  formulation 
can  then  be  easily  checked  against  how  it  changed  the 
net  result.  In  the  present  formulation,  we  have  de- 
termined, for  example,  that  the  inclusion  of  a complete 
droplet  distribution  function  compared  to  using  the 
mean  droplet  size  does  not  significantly  affect  the 
total  mass  accreted  for  most  problems  of  engineering 
interest  where  droplet  sizes  exceed  about  20-30/zm 
in  radius. 
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APPENDIX  A.  COMPUTER  PROGRAM 
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A facsimile  catalog  card  in  Library  of  Congress  MARC  format  is 
reproduced  below. 


Ackley,  S.F. 

Computer  modeling  of  atmospheric  ice  accretion  / by  S.F. 
Ackley  and  M.K.  Temp.-ton.  Hanover,  N.H.:  U.S.  Cold  Regions 
Research  and  Engineering  Laboratory;  Springfield,  Va.:  available 
from  National  Technical  Information  Service,  1979. 
iii,  39  p.,  illus.;  27  cm.  ( CRREL  Report  794.) 

Prepared  for  Directorate  of  Military  Programs,  Office,  Chief 
of  Engineers  by  Corps  of  Engineers,  U.S.  Army  Cold  Regions 
Research  and  Engineering  Laboratory. 
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